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Abstract 

The method of equivariant moving frames on multi-space is used to construct sym- 
metry preserving finite difference schemes of partial differential equations invariant 
under finite-dimensional symmetry groups. Invariant numerical schemes for a heat 
equation with a logarithmic source and the spherical Burgers equation are obtained. 
Numerical tests show how invariant schemes can be more accurate than standard dis- 
cretizations on uniform rectangular meshes. 

1 Introduction 

In modern numerical analysis, much effort has been invested into developing geometric 
integrators that incorporate geometrical structures of the system of differential equa- 
tions being approximated. Well known examples include symplectic integrators [BJ, 
energy preserving methods [23], Lie-Poisson preserving methods |27j . and symmetry 
preserving numerical schemes [H 1191 f26] . The motivation behind all this work is that, 
as a rule of thumb, geometric integrators give better results than many other standard 
numerical methods since they take into account qualitative properties of the system 
being studied. 

For ordinary differential equations, the problem of generating invariant numeri- 
cal schemes preserving the point symmetries of the original equations is now well- 
understood [2 HI [181 125] • There exists essentially two different methods for generating 
invariant finite difference schemes of differential equations. The first approach, mainly 
developed by Dorodnitsyn, Levi and Winternitz, is based on Lie's infinitesimal sym- 
metry method [31 [TU1 HH [T21 [T21 [2BJ . This approach makes use of Lie's infinitesimal 
symmetry criterion [20] to obtain finite difference invariants from which an invariant 
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scheme is constructed by finding a combination of the invariants that converges, in the 
continuous limit, to the original system of differential equations. The second approach, 
developed by Olver and Kim, consists of using the method of equivariant moving frames 
|16 | ll7 1 l2l]. With Lie's infinitesimal method, the construction of invariant schemes can 
sometimes require a lot of work and insights, on the other hand, with the moving frame 
method the construction is completely algorithmic. In both cases, the methods have 
proven their efficiency for generating invariant numerical schemes for ordinary differen- 
tial equations. In comparison, fewer applications involving partial differential equations 
can be found in the literature [3 EJ HDl EH [26] . 

In this paper we use the equivariant moving frame method to generate invariant 
finite difference schemes for partial differential equations with finite-dimensional sym- 
metry groups. To implement the moving frame method, the first step is to obtain 
appropriate approximations of the partial derivatives on an arbitrary mesh. In |23j, 
Olver proposes a new approach to the theory of interpolation of functions of several 
variables, based on non-commutative quasi-determinants, to obtain these approxima- 
tions. We show here that one can use standard Taylor polynomial approximations to 
achieve the same goal in a somewhat simpler way. The formulas obtained suggest an in- 
teresting interpretation of the continuous limit of the finite difference derivatives. More 
details are given in Section [3j 

Since the proper geometrical space for the numerical analysis of differential equa- 
tions is that of multi-space [21] , we begin with a review of multi-space in Section [2j 
The formulas used to approximate derivatives on an arbitrary mesh are introduced in 
Section [3} With this in hand, we review the multi- moving frame construction in Section 
[4] As with the standard moving frame method, the equivariant multi-frame construc- 
tion relies on Cartan's normalization of the group parameters of the action. Given a 
multi-frame, there is a canonical invariantization map which projects finite difference 
expressions onto their invariant finite difference counterpart. In particular, the invari- 
antization of the finite difference derivatives gives finite difference invariants which, in 
the continuous limit, converge to the normalized differential invariants [14j . Thus, by 
rewriting a system of partial differential equations in terms of the normalized differen- 
tial invariants, an invariant numerical scheme is obtained by replacing the normalized 
differential invariants by their invariant finite difference approximations. The method 
is illustrated in Section 5 where invariant numerical schemes are obtained for a heat 
equation with a logarithmic source and for the spherical Burgers' equation. Section 6 
is dedicated to numerical tests in which the precision of completely invariant schemes 
is compared to invariant and standard discretizations on a uniform mesh. 

2 Differential Equations and Numerical Schemes 

Let M be an m-dimensional manifold. For < n < 00 let J n = J n (M,p) denote 
the extended n th order jet space of 1 < p < m dimensional submanifolds S C M. 
The jet space is defined as the space of equivalence classes of submanifolds under the 
equivalence relation of n th contact at a point [20J. We use the notation ] n S\ z to denote 
the n-jet of S at z 6 S. Given a submanifold S C M, we introduce the local coordinates 
z = (2, u) = (x 1 , . . . , x p , u , . . . , u q ), with q = m — p, so that S is locally the graph of a 
function f(x): S = {(x, f(x))}. In this coordinate chart the coordinate jets of j n S are 
z (n) _ ( X)U ( n )) ; where denotes the collection of derivatives Uj = d k u a / (dx) J with 
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< k = #J < n. 

The approximation of the n-jet of a submanifold by finite difference derivatives 
is based on the evaluation of the submanifold at several points. For reasons which 
will become clearer in the next section, we label the sample points with a multi-index 
%N = (%n,un) = (%N > f{ x N)) where N = (n , . . . ,n p ) G W . For the finite difference 
derivatives to be well defined, the sample points must be distinct in the independent 
variables. Thus we introduce the fe-fold joint product 

M ok = , . . . jZNk )\x Ni ^ x Nj for all i ^ j} 

which is the subset of the A;- fold Cartesian product M xk . 

A smooth function A : J n — > M. on (an open subset of) the jet space is known as a 
differential function. A system of differential equations is defined by the vanishing of 
one or more differential functions. In the local coordinates z^ n ' = (x, system of 

differential equations is given by 

A x (x, u (n) ) = ■■■ = A e (x, «W) = 0. (2.1) 

Definition 2.1. Let A v (x, u^) = 0, v = 1, . . . , £, be a system of differential equations. 
A finite difference numerical scheme is a system of equations 

E fJ ,(zN 1 , . . -,z Nk ) = 0, (i = 1, ... ,£, ...,£ + I. 

defined on the joint product M ok with the property that in the coalescent limit (con- 
tinuous limit) zn x — > z, E^ZjVij • • • > z N k ) — > A„(x,«( n )) for /j, = 1, . . . ,£ and — > 
for fj, = £ + 1, . . . ,£ + 1. The last I equations = 0, . . . , Ei + i = impose constraints 
on the mesh. 

There are two restrictions on the equations Ei+i = 0, . . . , Ei + i = 0. Firstly, the 
equations must be compatible so that, provided appropriate initial conditions, the in- 
dependent variables xat. are uniquely defined. Secondly, these equations should not 
impose any restriction on the dependent variables ■ When these two conditions are 
satisfied we will say that the mesh equations are compatible. 



In Definition 2.1 the number of copies k in the joint product will depend on the 
order of approximation of the numerical scheme. The minimal number of sample points 
required to approximate the n th order system of equations (2.1) is k = <?( p ^ n ) = 



dim J n — p. More sample points can be added for better numerical results . 
Example 2.2. A standard numerical scheme for the heat equation 

A(x,i,u (2) ) = u t - u xx - ulnu = (2.2) 
with logarithmic source on the uniform rectangular mesh 

Xm,n = hm + xo, t m , n = kn + to, where h,k,xo,to are constants (2.3) 
and m, n are integers is given by 

El = A t u - A 2 x u - U m ,n In U m ,n = 0, (2.4a) 

where 

a U"m,n+1 ^"m,n A 2 'U j m+2,n 2u m _|_i n + Um,n 

A t u = , A x u = 



h 2 
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are the standard finite difference derivatives on the uniform rectangular mesh (2.3). In 



terms of Definition 2.1| the mesh (|2.3|) is defined by the equations 

-^2,3,4,5 



3? 771+1,71 



^777,77 



h, 



t 



m+l,7i 



tr. 



0, 



tm,n+l tm,n = k. 



(2.4b) 



Given the initial conditions xq, to the equations (2.4b) uniquely specify the points 

(^771,71) ^771, 7l) • 

Now, let G be an r-dimensional Lie group acting regularly on M. Throughout the 
paper, we will consistently use lower case letters to denote the source coordinates of 
the action, and capital letters to denote the target coordinates: 



X' 



9 ■ x , 



U a = g ■ u a , 



,P, 



a 



1.. 



geG. 



(2.5) 



Since a Lie group action preserves the contact equivalence of submanifolds, it induces 
an action on J n : 

9 ■ inS\ z = j n (g ■ S)\ g . z , g€G, (2.6) 
called the n th order prolonged action. In applications, the prolonged action (2.6) is 



obtained by implementing the chain rule. In local coordinates, we use the notation 
(X, U^) = g ■ (x, u^™)) to denote the prolonged action. The Lie group G also induces a 
natural action on the &:-fold Cartesian product M xk given by the product action G xk : 



g ■ (z Nl , . . . , z Nk ) = (g ■ z Nl , . . . , g ■ z Nk ). 



(2.7) 



Definition 2.3. Let G be a Lie group acting regularly on M ~ X x U and A u (x, u^) = 
a system of differential equations. The system of differential equations is said to be G- 
invariant if A u (g-(x,u^ n ')) = A u (x, u("') for all g S G. Similarly, let E^z^, zx k ) = 
be a system of finite difference equations. Then it is G-invariant if E^{g ■ (z^, . . . , 
z Nk )) = E^(z Nl ,. . . , z Nk ) for all g G G. 



Example 2.4. The heat equation with logarithmic source (2.2) is invariant under the 
group of transformations 



X = x + 2A 3 e* + A 2 



T = t + X 1 , InU = In u - X 3 e l x - \\e 2t + A 4 e 



2.i 



where Ai, . . . , A4 6 M, [20]. By direct computation, it is not difficult to see that the 
numerical scheme ( 2.4a| ), on the rectangular mesh (2.4b), is only invariant under the 
3-dimensional group of transformations 



X = x + X 2 , T = t + \ lt 
To see that the one-parameter group 

X = x + 2X 3 e\ T = t, 
is not a symmetry, it suffices to note that 



X 



m,?7+l 



x„ 



=2A 3 (e* m '" +1 
=2A 3 (e* m >" +1 



Ya.U = Inn + X^e . 



\nU = hiu- X^x- X\e 2t 



(^777,77 + 1 ^777,77) 

^ when A 3 / 0. 



As Example 2.4 shows, for a numerical scheme to preserve the symmetry group of 



a differential equation, a rectangular mesh might be too restrictive. Consequently, the 
theory of invariant numerical schemes must be developed over arbitrary meshes. 
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3 Finite Difference Derivatives 



In this section we obtain finite difference derivative expressions on an arbitrary mesh. 
There are different ways to do so. One could use the theory of multivariate interpolation 
[23], but we prefer to use Taylor polynomial approximations. Depending on the method 
used, the expressions for the finite difference derivatives can be slightly different, but 
this difference does not alter the implementation of the moving frame construction 
discussed in the next section. 

To simplify the notation, we assume that the generic point under consideration 
corresponds to the zero multi-index. Also, let 



a = (o,...,o,i,o,...,o), 



1, . . . ,p, 



be the multi-index of length p with 1 in the i th component, and zero elsewhere. Finally, 
let 

AiZ = Z 0+ei ~ Z = Z ez - Zq, 1 = 1, ...,p, 

be the usual forward difference operator in the i th component. 

To obtain finite difference approximations for the first order derivatives tt" we 
consider the first order Taylor polynomial expansions 



AiU, 



* X ' U X3> 



,P, 



a 



,q- 



(3.1) 



Solving for the first order partial derivatives u a t we obtain 



a,d\ 



a,d 

xP 



/Ai: 



yAp^o 



A r p 
P / 



r 









a 



l,...,q, (3.2) 



Ax- 1 



which are well defined provided the matrix Axo is invertible. When this is the case, 
the expressions (3.2) define first order finite difference derivatives on an arbitrary mesh. 
To obtain finite difference approximations for the second order derivatives, it suffices 
to consider the second order Taylor polynomial expansions 



AjAjiio « £ AjAjXg -< fe +E [(4 

k,l=l 

AjXg • AjXq] — 



x o)( x ei 



+e 7 - 



4) 



k=l 



(3.3) 



with 1 < j < i < p and a = l,...,q. Substituting the first order finite difference 
derivatives (3.2) in (3.3) and solving for the second order derivatives u a kxl gives finite 
difference approximations for second order derivatives. Those are well defined expres- 
sions provided that the matrix with coefficients given by the factors in front of the 
second order derivatives in (3.3) is invertible. Continuing in this fashion, it is possible 
to obtain finite difference approximations for third order derivatives and so on. 

We now consider the continuous limit of those finite difference derivative expressions 
and verify that they converge to the corresponding partial derivatives. The standard 
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way of taking the continuous limit of finite difference quantities is to assume that all 
nodes z^ i coalesce to the same point z, and the role of the multi- index Ni is simply to 
label the points. In the following, we give a more important role to the multi-index iVj. 



Instead of viewing the solution of a system of partial differential equations (2.1 ) as the 



graph of a function (x,u(x)), we can view a solution as a parametrized p-dimensional 
submanifold 



S->M~Ix[/, 



(s l ,...,sn^z(s) = (x(s),u(s)) 



transversed to the fibers {c} x U. In this setting, the distinction between the inde- 
pendent variables x % and the dependent variables u a disappears. Moving to the finite 
difference picture, we now view the multi-indices Ni as forming a square grid with edges 
of length 1 in the parameter space S. Then, a point Zfq i £ M is just the image of z(s) 
at s = Ni, see Figuref [T] (in two independent variables). 



(0,1; 



(0,0) 



XD 



Zq,1 



^0,0 



;i,o) 



Zl,0 



Figure 1: Multi- index interpretation. 



In the continuous limit, the point = z(N{) converges to zq = z(0) = z by 
coalescing the multi-index Ni to the origin. With this point of view the differences 
z ei — zq converge to z s % since 



\Z = Z et - Z 



zq 



a' 



(T l -S>0 



i = l,..., p. 



Similarly, in the continuous limit 

AjAj2o > z s i sj \ z , i,j = l,...,p, 



and so on. It thus follows that the expressions (3.1) converge to 

p 



a = l,...,q, i = l,...,p, 



which is the chain rule formula for u = u(x(s)). Similarly, the expressions (3.3) converge 
to 



(3.4) 



k,l=l 



in the continuous limit. 
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Example 3.1. The above discussion is now specialized to the particular situation of 
two independent variables (x,y) and one dependent variable u = u(x,y). Without loss 
of generality, all expressions are centered around zo,o = (xo,o> 2/0,0; ^o,o)> an d neighboring 
points are denoted by 

Zm,n — \P^m,ni Vm,m ^m,n)> 771,71 £ Z. 

Also, let 

denote the standard forward difference operators in the two indices and let (s^s 2 ) = 
(s, i) be coordinates for the parameter space S 1 . 

The indices involved in the expressions of the first and second order discrete partial 
derivatives are displayed in Figure 2. In general, the definition of the n th order discrete 
derivatives will involve the indices contained in the right triangle formed by the origin 
and the vertices (n, 0), (0,n). 



it t 




(a) First order derivatives (b) Second order deriva- 

tives 



Figure 2: Indices used to define first and second order discrete partial derivatives. 



Following our general procedure, to obtain the first order finite difference approxi- 
mations of u x and u y on an arbitrary mesh, the two Taylor expansions 

u s « A-u 0j0 ~ Axo.o • Ux + Ay ,o • u y , u t « 5u ,o ~ Sx 0fi ■ u x + <5y ,o ■ u y 
are considered. Solving for u x and u y we obtain 



5y ,o • An ,o - Ay ,o • Su . 



o 



Ax ,o • Sy fi - 5x 0fi • Ay . 



o 



Uy ~ Uy 



Axpfl ■ Supfl - 5xpfl • An 0i0 

Ax ,o • <5y ,o - ^o,o • Ay ,o ' 

(3.6) 



In the continuous limit the expressions (3.6) converge to 

y t u s - y s u t 



u. 



x s u t - x t u s 



x s yt ~ x t y s x s y t - x t y s 

which are the usual formulas obtained by the chain rule if x = x(s,t), y = y(s,t). 



(3.7) 
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To obtain approximations for the second order derivatives, it suffices to take second 
order Taylor expansions: 



Us 



A 140,0 ~ A x ,o • u x + A 1/0,0 • u y + [(x 2 ,o - xq,o) - 2(Ax ,o) ]- 
+ [(^2,0 - £o,o) (2/2,0 - 2/o,o) - 2Ax ,o ■ A2/o,o]«aiy 
+ [(y 2 ,o-yo,o) 2 -2(A yo ,o) 2 ]^, 



u xx 



u st ~ 5Au , » <5Ax ,o • i*x + £Ay ,o • u y + - x ,o) - (Ax ,o) - (&co,o) 

+ - x ,o) (2/1,1 - 2/o,o ) - Ax ,o • Ay ,o - Sx ,o ■ Syo fi }u xy (3.8) 
+ [(yi,i-yo,o) 2 -(Ay ,o) 2 -(^o,o) 2 ]^, 



uu « 5 2 u ,o ~ <5 2 ^o,o • «x + £ 2 yo,o • ""y + [(^0,2 - ^o,o) 2 - 2(5x ,o) 2 ]^ 
+ [(^0,2 - £0,0X2/0,2 - 2/o,o) - 2<5x ,o ■ <fyo,o]uxj/ 
+ [(yo,2-2/o,o) 2 -2(5yo,o) 2 ]^. 



Solving for the second order derivatives u xx , u xy , u yy in (3.8), and replacing the first 
order derivatives u x u y with the approximations (3.6), the expressions for the discrete 
second order derivatives are 

(u xx \ ( u xx \ 
u*v ~ <y = H ~ ly i (3.9) 

where V is the column vector 

/ A 2 u , - A 2 x ,o • u£ - A 2 y ,o ■ 
V = 5Au 0fi - <5Ax ,o • u d x - 5Ay 0fi ■ u\ 
\ 5 2 u 0fi - 5 2 x ,o ■ u d x - 5 2 y 0fl ■ u d y 

and H is the 3x3 matrix with entries 

#ii = [(^2,0 - ^o,o) 2 - 2(Az ,o) 2 ]/2, #1,2 = (ac2,o - ^0,0X2/2,0 - 2/o,o) - 2Ax ,o ■ Ay ,o, 
#1,3 = [(2/2,0 - 2/o,o) 2 - 2(Ay ,o) 2 ]/2, #2,1 = [(a>i,i - x ,o) 2 - (Ax , ) 2 - (<5x ,o) 2 ]/2, 
#3,3 = [(2/0,2 - 2/o,o) 2 - 2(,5yo,o) 2 ]/2, #2,3 = [(2/1,1 - 2/o,o) 2 - (Ay , ) 2 - (<%o,o) 2 ]/2, 
#3,1 = [(^0,2 - ^o,o) 2 - 2(fe , ) 2 ]/2, #3,2 = (x ,2 - £0,0X2/0,2 - 2/o,o) - 2(5xo,o ■ #2/0,0, 
#2,2 = (xi,i - £0,0X2/1,1 - 2/o,o) - Ax ,o • Ayo,o - &c ,o • <fyo,o- 

Continuing in this fashion, it is possible to obtain higher order finite difference deriva- 
tives on an arbitrary mesh. 

The finite difference derivatives constructed above are first order approximations 
of the continuous derivatives. More accurate approximations can be obtained by using 
more nodes to approximate the derivatives. Also, those approximations are constructed 
solely with the forward difference operators Aj. To obtain centered difference approxi- 
mations one can use the right and left differences operators 

For better numerical schemes, centered difference approximations are used in the next 
sections. 
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4 Equivariant Moving Frames 



In this section we review the moving frame construction on joint space. The exposition 
follows [21]. 

Definition 4.1. Let G be a finite-dimensional Lie group acting smoothly on a manifold 
M. A right moving frame is a smooth G-equivariant map p: M — > G such that 

p(g ■ z) = p(z) ■ g~ x for all z € M, g £ G. 

Theorem 4.2. A right moving frame exists in a neighborhood of a point z £ M if and 
only if G acts freely and regularly near z. 

The action is free if at every point z G M the isotropy subgroup G z = {e} is trivial. 
Under this assumption, the group orbits are of dimension r = dim G. The action is 
regular if the orbits form a regular foliation. The moving frame construction follows 
from the following theorem. 

Theorem 4.3. If G acts freely and regularly on M, and K, C M is a cross-section to 
the group orbits, then the right moving frame p: M — > G at z £ M is defined as the 
unique group element g = p{z) which sends z to the cross-section p(z) ■ z G /C. 

While it is not necessary, we assume K. = {z 1 = c , . . . , z r = c r } is a coordinate 
cross-section obtained by fixing the first r coordinates z = (z , . . . , z m ) of M to some 
suitable constants. The moving frame p(z) is then obtained by solving the normalization 
equations 

Z x = g-z 1 = c 1 , ... Z r = g-z r = c r , (4.1) 

for the group parameters g = (g±, ... , g T ) in terms of z = (z , . . . , z m ). Given a moving 
frame, there is a systematic way of associating an invariant to a function. 

Definition 4.4. Let pbea right moving frame, the invariantization of a scalar function 
F : M — > M is the invariant function 

I(z) = i(F)(z) = F(p(z).z). (4.2) 

Proposition 4.5. The invariantization of the coordinate functions i(z r+1 ), . . . , i(z m ) 
provide a complete set of m — r functionally independent invariants on M. 

Note that by the moving frame construction, the invariantization of the coordinates 



defining the normalization equations (4.1) are constant, l{z 1 ) = c , . . ., t(z r ) 



Also, if I(z) is an invariant then = I. 

For most groups of interest, the action on M fails to be free. There are two common 
methods for making an (effective) group action free. In geometry, this is accomplished 
by prolonging the action to a jet space J n of suitably high order. The cross-section 
fC is then an r-dimensional submanifold of J n and the normalization of the group 
parameters yields a standard moving frame p^ : J n — > G. The invariantization of the 
jet coordinates l(z^) = 1^ gives a complete set of differential invariants [T3]. The 
second possibility is to consider the product action of G on a suitably large joint product 
M ofc C M xk . The moving frame construction yields the product frame p ok : M ok — > 
G, and the invariantization of the coordinate functions z^ i gives a complete set of 
functionally independent joint invariants |22j. 
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In |21j , it is shown that the right geometrical space to study symmetry of numerical 
schemes is that of multi-space. The moving frame construction yields what is called a 



multi-frame and the invariantization map (4.2) gives multi-invariants. For the present 



discussion, it is enough to know that in the continuous limit p 0k — > p^ provided the 
moving frames are compatible. To obtain compatible moving frames, instead of working 
with the standard product action Zx l = (Xx t , U^) = g-z^ l on M ofc we should consider 
the joint action on the the finite difference coordinates 

x ) X N 1 > • • • ) x N k i u i u ) 

where denotes the collection of finite difference derivatives u a j d of order < n on 

an arbitrary mesh computed in Section [3J For the discrete invariants (Hq, I d '^) = 
i(xQ,u d '^) to converge to the differential invariants (H,I^) = i(x,u( n >), the cross- 
section lC k defining the product frame must, in the continuous limit, converge to a 
cross-section JC^ for the prolonged action on J n . Assuming the action is transitive on 
the independent variables x, this constraint on )C ok implies that we can only normalize 
the independent variables at one multi-index, say Xq = c, |21j . 



Example 4.6. To illustrate the above discussion we consider the symmetry group of 
the heat equation with logarithmic source: ut = u xx + ulnu. The prolonged action 



induced by (2.8) is obtained by implementing the chain rule. The result is 

(lnU) T = (lnu) 4 -A 3 e*x + A4e'-2A 3 e*(lnu) :r (In U) x = {\nu) x - A 3 e*, 

(hxU)xx = (hiu) xx , 

and so on. Choosing the cross-section 

/CW = { x = t = Inu = (lnu) x = 0} 

and solving the normalization equations X = T = In U = (In U)x = for the group 
parameters Ai, A2, A3, A4 yields the moving frame 



Ai = -t, A 2 = -{x + 2(\nu) x ), X 3 = e t (\nu) x , 
A4 = e~*(— Inn + x(lnu) x + Qnu) x ). 



(4.4) 



By construction 

l(x) = t(t) = l(u) = i(ln(n) x ) = 0, 
and the invariantization process yields the differential invariants 

I = i((ln u) t ) = In u t - Inu- (In u)l, J = l((Iiiu) xx ) = (lnu) xx , (4.5) 

and more. 

We now repeat the computations for the discrete case. In the following, we assume 
that 

h,o ~ to,o = (4.6) 



as it simplifies the calculations. The equality (4.6) is compatible with the group action 



(2.8) as it is an invariant equation of the joint action. To obtain a product frame 



compatible with (4.4) we choose the cross-section 

1C 1 = {x ,o = io,o = ^0,0 = (lnn)f = 0}, 
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where, due to the assumption (|4.6|), (lnu) 



In ui o — lnu_i 



-. Solving the corre- 



21,0 - £-1,0 

sponding normalization equations for the group parameters Ai, A2, A3, A4 gives the 
product frame 



Ai — —toft, 
A 4 = 



A 5 



-to,o 



-(z ,o + 2(lnu)!*), X 3 = e- t ^(\nu) d x , 



d\2 



ln« ,o + x 0)0 (lnu) x + ((lnit)£) 



(4.7) 



The invariantization map is then completely defined and yields, among many others, 
the discrete invariants 



r 



t((lnn 



lnupi - e T lntt 0i0 _ cr T 
r 



■It 



T 

J d = L((\nu) d xx ) = Onu) d xx , 



e T (\nu) d x + - e —({\nu) d x f, 



where a = xq 1 — xq and r = to 1 — *o 0, and 



(lnn)l 



x\,o - a; -i,o 
The discrete invariant 



lnm ;0 - ln-u 0i0 
£1,0 - ^0,0 



lntt 0i0 - hm-1,0 

XOfi ~ £-1,0 



(4.8) 



i(x ,i) = a + 2(e T - l)(lnu)£ 



(4.9) 



will also prove to be useful in the construction of an invariant numerical scheme in the 
following section. 



5 Invariant Numerical Schemes 

Let A u (x,u^) = be a G-invariant system of differential equations as defined in 
Definition |2.3[ The invariance of the system guarantees that it can be expressed in 
terms of the normalized invariants l(x,u^) = (H,I^). Using the invariantization 



map (4.2) 



= A u (x,u^) = l(A u (x,u^)) = A V {H, JW). 



(5.1) 



To obtain an invariant numerical scheme we simply need to replace the differential 



invariants (H,I^) in (5.1) by their invariant discrete counterparts l( 



x ,W 



d,(n) 



(Hq, 7 rf '( n )) and add equations specifying the mesh. If the equations determining the 
mesh are invariant under the group action G, the scheme is said to be fully invariant, 
otherwise it is called partially invariant. In applications, invariant constraints on the 
mesh are usually specified with the aim of simplifying the expressions of the invariant 



scheme. This is illustrated in Examples |5.1| and 5.2 



An alternative way of constructing the same invariant numerical scheme as above 
consists of replacing the partial derivatives in the differential equations A u (x, u^) = 
by their finite difference approximation 



A u {x,u (n) ) = 



A u (x 0fi ,u d ^) = 0, 



(5.2) 



followed by the invariantization of (5.2) 
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Example 5.1. The heat equation with source (2.2) is easily expressed in terms of the 
differential invariants (4.5). By invariantizing (2.2), the equation can be rewritten as 

i(ut — u xx — u In u) = I — J. 







(5.3) 



To obtain an invariant discrete version, I and J in (5.3) are simply replaced by their 
discrete counterparts: 



jd_ja 



In uq i — e T In uq q o 



'It 



((ln^) 2 -(ln^ = 0. (5.4) 



By construction, the numerical scheme is a valid approximation of (2.2) on any mesh 
satisfying the flat time assumption (4.6). At the moment, there is no systematic method 
for obtaining invariant equations defining the mesh of an invariant numerical scheme. 
In |17j it is claimed that the equations for the invariant mesh are simply obtained 
by invariantizating the equations of the mesh equations for the non-invariant numerical 
scheme. The issue with this proposition is that there are no guarantee that the equations 
obtained are compatible. Indeed, this idea was used in [7] and led to the incompatible 
mesh equations (83-85). To illustrate the problem we invariantize the equations (2.4b) 
describing a uniform rectangular mesh. The result is 



£i,o - ^0,0 = h, 
h.o — too = 0, 



%0,1 — ^0,0 = cr = 2(lnn)^(l — e T 
*o,i - *0,0 = T = k. 



In the time variable t the mesh equations are compatible but this is not the case in the 
spatial variable x. By shifting the first equation by S and the second equation by A 
(recall (3.5)) we obtain the constraint 

h = xi,i - x ,i = 2A[(lnu) x ](l - e k ). (5.5) 

Since h and k are constants, the equation ( |5.5| ) imposes a restriction on the solution 
u which is not admissible. To circumvent the problem we can neglect the equation 
x\o — ^0,0 = h and assume that the mesh equations are given by 

a = 2(1- e T )(lnu)i, t = k, At = 0. (5.6) 



With (5.6) the mesh is uniquely determined once the sample points in x are fixed at 
some time t = to- In conclusion, a fully invariant numerical numerical scheme for the 
heat equation (2.2) is given by 

e T In u ,o T 1 



Ei 



In no 



-((ln^) 2 -(lnn)^ = 0, 
= r - k = 0, E 4 = At 



E 2 = a - 2(1 - e T )(lnu)i = 0, 



T 

E 3 



(5.7) 



0. 



Example 5.2. In this second example we find an invariant numerical scheme for the 
spherical Burgers equation 



u 

u t + - + uu x + u xx = 0, 



t > 0. 



(5.t 



The equation (5.8) admits the three-dimensional Lie algebra of infinitesimal symmetry 





generators, [To] . 



vi 



d_ 

dx ' 



v 2 



— 2t— 
dx dt 



d _ d Id 

du' 3 dx t du 
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The corresponding group action is 
X = e A2 (x + A 3 lnt) + Ai, 



T = e 2X H 



U = e 



u + 



(5.9) 



For the joint action, a simple choice of cross-section is 

£0,0 = 0, t ,o = 1, u 0fi = 0. 

Solving the normalization equations 

Ao,o = 0, Tqa = 1, Uofi = 

for the group parameters yields the product frame 



Ai 



^0,0 ~ "0,0*0,0 hi*o,i 



A 2 = In 



1 



A3 = -1*0,0*0,0- (5.10) 



An invariant numerical scheme for the spherical Burgers equation (5.8) is given by 
the invariantization of the discrete derivatives u xx and u d : 



= I d = i{ui) + i(u d xx ) = u d + -^ + u 0fi u d x + u d xx . 

to,o 



(5.11) 



The expression ( 5.11 ) is defined on any compatible mesh. We now impose some invariant 
constraint on the mesh which will simplify the coordinate expressions of (5.11). Once 
more, we assume that Ato,o = *i,o — *o,o = as this constraint is invariant under the 
group action (5.9 ). Another invariant constraint is given by <5 2 to,o = *o,2 _ 2io,i+*o,o = 0. 
Finally, using the invariant 



t(«o,i) = £0,1 - £0,0 - tit In 



*o,i 

*0,0 



a ~ ""0,0*0,0 In 



*0,1 
*0,0 



an invariant mesh is specified by the equations 



o = ""0,0*0,0 In 



*0,1 
*0,0 



At ,0 = 0, <5 2 t ,o = 0. 



(5.12) 



As in the previous example, once the sample points in x are determined at some time 
t = to and the step size in r = k is chosen, the equations (5.12) uniquely fix the mesh. 
On the mesh (5.12) the equation (5.11 ) simplifies to the fully invariant numerical scheme 



T-i ""0,0*0,0 . „ d n 
£-1 = «0 1 1" 2<TU XX = 0, 



E 2 = a - u ,o*o,o In 



*o,i 

*0,0 



*0,1 



0, E 3 = At n = 0, E A = 6 2 t 0fi = 0, 



(5.13) 



where 



u , 



#1,0 - 2-1,0 



"1,0 - ^0,0 
Xlfl - x ,o 



^0,0 - "-1,0 
XQfl - 2-1,0 
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6 Numerical Tests 



Fully invariant numerical schemes for the heat equation with logarithmic source (2.2) 



and the spherical Burgers equation (5.8) were obtained in (5.7) and (5.13), respectively. 
To illustrate how preservation of symmetry can increase the accuracy of a numerical 
method, we compare the fully invariant schemes with two closely related schemes. 



Firstly, recall that provided the flat time assumption (4.6) holds, the discretizations 



(5.4) and (5.11) are valid approximations of the corresponding partial differential equa- 
tion on any compatible mesh. To gauge the effect of imposing invariant constraints 
on the mesh, we compare the fully invariant schemes to partially invariant schemes 



obtained by restricting (5.4) and (5.11) to non-invariant uniform rectangular meshes. 
Then, for further comparison, the fully and partially invariant schemes are set against 
standard finite difference approximations on uniform rectangular meshes. 



6.1 Heat Equation 

For the heat equation (|2.2l) we used the steady-state solution 



u(x, t) = exp 



ce 



x z 1 
T + 2 



with 



0, 



(6.1) 



to compare the precision of the three schemes. Starting at t = and working on the 
space interval [—5,5], the solution after one unit of time elapsed is shown in Figure [3| 
The numerical simulation was done using an initial spacial step size of h = 0.15 and a 
constant time step r = k = 0.001. We note that for the fully invariant scheme (5.7) 
the mesh evolves according to the equation 

The result is an expansion of the space interval [—5, 5] over the course of the simulation 
as shown in Figure [4j The absolute errors for the three numerical schemes are given 



1=1 Fully invariant 
Partially invariant 
* Standard 
Exact Solution 
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gure 3: Solution u = exp[a; 2 /4 + 1/2]. 



Figure 4: Invariant mesh. 



in Figure [5} As anticipated, the standard numerical scheme offers the worst precision 
among the three. The partially invariant scheme on a rectangular mesh gives better 
results, but the error decreases even more for the fully invariant numerical scheme. 
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Fully invariant 
Partially invariant 
Standard 




Figure 5: Absolute errors for the solution u = exp[x 2 /4 + 1/2]. 



6.2 Burgers Equation 



For the spherical Burgers equation (5.8), a numerical test was performed using the 
solution 

x ~\~ C\ 

u(x,t) = — — - with c\ = and C2 = 1. 

t{C2 + mt) 

At t = 1, the initial step sizes chosen are h = 0.5 in x and r = k = 0.001 in t. In Figure 
[6j the solution is shown at t = 1.5. As in the previous simulation, the mesh of the fully 
invariant numerical scheme evolves as a function of time. The evolution is governed by 
the equation 

^0,1 



& = "0,0*0,0 In 



*0,0 



The evolutive mesh appears in Figure [7] 

The absolute errors for the three numerical schemes are given in Figures [8] and [9} 
Once more, the fully invariant scheme is more precise. For the solution considered the 
improvement is considerable as the error is reduced by a factor 10~ 13 compared with the 
two other schemes. On the other hand, while the partially invariant scheme is slightly 
more accurate than the standard scheme, the errors are comparable. 



7 Concluding Remarks 

Other simulations, not included in the paper, indicate that invariant schemes do not 
always improve the numerical accuracy of a standard scheme. While the error is gener- 
ally not worst than the standard scheme it is still not clear for which types of equations 
or solutions an invariant scheme will give significantly better results. Knowing when 
this is the case would be useful for future applications. 

As with all other works on the subject, the important question of the stability of an 
invariant scheme was not addressed in this paper. At the moment it is not clear how 
the evolution of the mesh affects the stability. 

Finally, many partial differential equations admit infinite-dimensional symmetry 
groups. Classical examples include the Kadomtsev-Petviashvili equation [8], the Infeld- 
Rowlands equation [13] and the Davey-Stewartson equations [5]. An interesting prob- 
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Figure 6: Solution u = x/[t(l + licit)}. 



Figure 7: Invariant mesh. 
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Figure 8: Absolute error for the fully invariant 
scheme. 



Figure 9: Absolute errors for the standard and 
partially invariant schemes. 



lem consists of finding a procedure for obtaining invariant discretizations of partial 
differential equations admitting infinite-dimensional symmetry groups. 
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